Beyond standard two-mode dynamics in Bosonic Josephson junctions 
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We examine the dynamics of a Bose-Einstein condensate in a symmetric double- well potential for 
a broad range of non-linear couplings. We demonstrate the existence of a region, beyond those of 
Josephson oscillations and self-trapping, which involves the dynamical excitation of the third mode 
of the double-well potential. We develop a simple semiclassical model for the coupling between 
the second and third modes that describes very satisfactorily the full time-dependent dynamics. 
Experimental conditions are proposed to probe this phenomenon. 
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I. INTRODUCTION 

Clouds of cold bosonic atoms exhibit a variety of quan- 
tal effects on the mesoscopic scale. The initial devel- 
opments dealt with dilute weakly interacting gases, and 
most of the experimental results could be explained by 
means of the Gross-Pitaevskii (GP) equation [H-Q. In 
recent years, the possibility of controlling the atom-atom 
scattering length through Feshbach resonances opened 
a broad range of possibilities for cold atomic clouds, 
connecting the physics of cold atoms to the physics of 
strongly interacting systems; for recent reviews see 

Here we consider Bosonic Josephson junctions (BJJ) 
as our benchmark. The GP approach successfully pre- 
dicted the existence of Josephson tunneling phenomena 
in clouds of bosonic atoms confined in a double- well po- 
tential @,[3. The physics of BJJ can be well understood 
by further assuming a two-mode description of the con- 
densate wave function, a simplification which correctly 
explains the existence of self-trapped states in BJJ when 
the non-linear interaction strength is increased J7j . This 
prediction was later confirmed experimentally [9]. Very 
recently, most of the semiclassical predictions of this two- 
mode approach, dealing with the Rabi to Josephson tran- 
sition, have been confirmed in an internal Josephson ex- 
periment [To( | . 

However, it is worth noting that the regime of appli- 
cability of the quantized two mode approximation can 
extend further: Recent examples are the experiments on 
BJJ, the production of numbe r sq ueezed states, and a 
non-linear atom interferometer [ill. Il2j. These phenom- 
ena are beyond GP, as they involve entangled states of 
the atoms in the cloud, but can, however, be explained 
by a requantization of the two-mode equations of GP, the 
Bose-Hubbard model [l3j |. Thus, they are still two-mode 
physics, albeit rcquantized. 

We scrutinize here the predictions of the GP equation 
when the non-linear interaction term is increased to set 
the condensate beyond the limits of validity of the usual 
two-mode truncation. Our focus is on a specific dynami- 
cal configuration, the evolution of the BEC when initially 
the majority of the atoms are located in one of the wells. 



By studying the oscillations of the population imbalance 
we find that increasing the non-linear coupling the am- 
plitude starts to increase, departing from the usual self- 
trapping behavior. We demonstrate that this dynamics 
can still be explained in terms of a few modes of the ef- 
fective GP potential, but that now it is the variation of 
coupling between the second and third modes that suc- 
cessfully explains the GP results. 

The manuscript is organized as follows, first in sec- 
tion |TT] the time dependent GP equation and the usual 
two-mode equations are presented, in section Mil we 
present the results obtained increasing the number of 
atoms, comparing the full time-dependent solutions of 
the GP equations with two-mode predictions. Section HVl 
contains a summary and conclusions. 



II. THEORETICAL DESCRIPTION 

For definitcness, we consider a dilute gas of N atoms at 
zero temperature and perform our study in ID. Thus our 
results should be relevant to cigar-shaped quasi-lD ex- 
periments. Setting H = m = 1 the evolution is described 
by the GP equation, 



1 d 2 
2'dx 2 



^(x,t)+V oB [^(x,t)]ip(x,t). (1) 



tj}{x,t) is normalized to 1, and the time dependent effec- 
tive potential is V e $[i/j(x, t)\ = V(x) + X N\i()(x, t)\ 2 . The 
relevant parameter in the GP equation is gw = A 7V, 
which sets the importance of the non-linear term. In 
this analysis we consider a repulsive interaction, Ao > 0. 
Different values of N and Ao produce exactly the same 
GP evolution provided AniV is kept fixed. The confin- 
ing double-well potential, V(x), is generated by V±(x) = 
(x ± 2) 2 connected by V p (x) = 3(1 — x 2 ) in the interval 
|x|<0.5. 

Eq. flT} is expected to provide accurate results for large 
enough number of atoms. The recent calculations re- 
ported in Ref. [l4| show how the exact ID dynamics of 
the BJJ approaches the GP predictions as the number of 
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FIG. 1: Maximum amplitude of the imbalance oscillations, 
•Zmax — Zmi n computed with GP (circles) as a function of gw- 
The solid (red) line is the classical two-mode prediction using 
Eqs. J2J. The dashed (blue) line is a two-mode calculation 
using modes (2,3) as explained in the text. 




FIG. 2: (color online) Effective potential V e g(x,t) = V(x) + 
gm\ i S(x, t)\ 2 for different values of giD- The bands are gen- 
erated by joining V o ff(a;,0) and V e s(x, t Zmia ), as explained in 
the text. The key to the various lines is shown in each panel. 



atoms is increased l . 

The standard two-mode approximation 0, HI, [H[ 
is obtained by making the following ansatz for the 
wave function: ^(x,t) = $ L (x) y/NiJ£)e l ' t ' L W + 
$fl(x)v / ^(i)e^- R W. The left and right modes ($ L (R) = 
l/\/2($i ± $2)) are constructed from the ground state, 
<I>i, and the first excited state, $2, of the stationary GP 
equation. These in turn are obtained numerically by an 
imaginary time evolution method. Then a set of equa- 
tions relating the population imbalance, z(t) = (Ni,(t) — 
Nn(t))/N and the phase difference, 8<f>(t) = (f>n(t) — 4>L(t) 
are easily derived 2 , 



m 

Hit) 



-2/C-v/l - z 2 {t) sm 5<j){t) 
z(t) 



(2) 



NUz(t) + 2JC 



cos S(j>(t) . 



where 



U = A dx<$> 



K = - dx((l/2)d x $ L d x $ R + $ L V(x)$ R ) (3) 



and integrals involving mixed products of $>l&r of order 
larger than one are neglected. Some authors characterize 
the dynamics with the variable A = NU/(2fC). In our 
double- well we find, using the modes computed at gw = 
0, K = 7.9 x 1(T 3 and U/X = 0.47, which gives A = 
29.7Arj-/V. The critical value of A c to have self-trapping 
within this two mode approximation Q, for (z(0) = 1, 
6(f>(0) — 0), is A c = 2, which translates into a critical 



1 The authors of Ref. Il4l consider up to 100 atoms, far away 
from the number of atoms considered in this work and in the 
experimental set up of Ref. [ij, N > 1000. 

2 These are the so-called standard two mode equations. Similar 
conclusions are obtained by using the improved two-mode equa- 
tions of Ref. Hal for the same double-well. 



value for g[^ = 0.067, where the superscript (1, 2) refers 
to the states involved in the tunneling dynamics. 



III. RESULTS 



In all the calculations, except those in section IHI CI 
we have if)(x,t = 0) = &l(x) (obtained by previously 
computing $1 and $2 of the GP for each gw), corre- 
sponding to the case of all atoms being on the left well, 
2(0) = 1, at t = 0. We study the dynamics for increasing 
values of gw > 0, going from the Rabi regime, gw = 0, 
to the Josephson and self-trapped dynamics, and further 
beyond the range of validity of the usual two-mode ap- 
proximation. 

Results are shown in figure [TJ There we compare the 
numerically determined GP amplitudes of z(t), empty 
circles, with the semiclassical two mode prediction of 
Eq. ([2]), solid (red). At gw < 9w there is no self- 
trapping and z(t) oscillates between +1 and —1, thus 
leading to a constant maximal amplitude, A z = z max — 

Zmin = 2. With increasing interaction strength, near 

(1 2) 

g\jj , self-trapping appears, the atoms become increas- 
ingly confined in the left well, and A z decreases abruptly. 
In all this range the semiclassical model predictions are 
very successful, covering the well known Josephson and 
self-trapped regimes. This range of gw is the one re- 
cently explored experimentally in Ref. (loj . 

With further increase of gw, deviations begin to ap- 
pear. Whereas the semiclassical two-mode model pre- 
dicts a smooth decrease of A Zl the GP calculations, 
empty circles, show a smooth reappearance of tunnel- 
ing between the two wells. This is the new phenomenon 
that we will now discuss and interpret. 

Figure [5] shows the effective potential defined in Eq. ([T]) 
for several values of gw at two different times, t = and 
t = £ 2min which correspond to the time of the first mini- 
mum of the population imbalance. In this way the band 
covers the variation of the effective potential during the 
simulation. When gw fS 9w the non-linear contribu- 
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FIG. 3: (a) First four eigenvalues of V e g at t = and t = t Zxnin , thus generating the bands. The dashed black lines are 
the first order perturbation theory calculation described in the text. The dotted line is the average energy of the initial state, 
{ip(x, 0)\H\ip(x, 0)). (b) Energy differences E2 —E\ (violet band) and £3 — E2 (green band) compared to twice the frequencies of 
oscillation of the population imbalance computed with GP, empty circles. The solid (red) line is the usual two-mode calculation 
using Eqs. The dashed (blue) line is a two-mode calculation using modes (2,3) as explained in the text, (c) and (d) depict 
the two relevant eigenstates of the effective potential (computed at t = 0) entering into the dynamics for gm =0.1 and 3, 
respectively, (e), (f), and (g) depict the normalized Fourier transform of z(t), \z(u)\, for g\o = 0.4, 1.4 and 3 respectively. 



tion is fairly small, and V e g(x,t) ~ V(x) at all times. 
This corresponds to the Rabi and Josephson regimes, 

with maximal oscillations of the population. Increasing 

( 1 2) 

gin further, g\jj < giu < 3, the value of V c g(x,t) in 
the left well is increased, but leaving the value in the 
right well almost unchanged. This is a direct consequence 



A. Analysis of the effective potential 

To clarify the discussion and further gauge the small- 
ness of the changes in the potential with time, we have 
determined the first four stationary solutions of the 
Schrodinger equation built with V Q g(x, t), with t a fixed 
parameter. Figure [3] (a), shows the eigenvalues thus 
found numerically. Two values of the parameter t have 
been chosen: t = and t — i Zmin , to form a band for 
each eigenvalue. The band width is seen to remain very 
small for gio < 2.2, in line with the largely time inde- 
pendent V e ff depicted in Fig. [5] (a). For small interaction 
strengths, gn> < 2, two of the eigenvalues are practi- 
cally independent of giD- they correspond to states 1 



of self-trapping. In this regime, the effective potential 
changes very little with time, see fig. [2] (a). Further in- 
creasing giD, 3 < giD ^5 5, the potential on the left well 
increases, and V e s does begin to change appreciably with 
time. Still, the dynamics remains self-trapped but with 
larger oscillation amplitudes in z(t). 



I 

and 3, mostly located in the right well, which remains 
unchanged as seen in Fig. [5] (a). Their values are thus 
close to the corresponding non-interacting harmonic os- 
cillator, E\° = V2/2 and E%° = 3^2/2. 

The other two eigenvalues increase smoothly with the 
interaction strength, they correspond to states 2 and 4 
located mostly in the left well. The increase follows the 
behavior of the effective potential shown in figure [2] The 
almost linear increase in energy of these two eigenval- 
ues can be understood treating the non-linear term as a 
perturbation, then 

6E2 = giD(rt°\\i>i°\ 2 = 9id/(2 1/4 V^) , 
SE 4 = g 1D (^°\ |VM 2 1^2°) = SE 2 /2 , (4) 



4 



where ip^ ° 2 are the first two eigenfunctions of the har- 
monic trap, V(x) — (x + 2) 2 : 



i>%°(x) = 2^\x + 2)^°{x). (5) 

These estimates are shown in Fig. |3j (a) as dashed lines. 

As sketched above, in the conventional two mode 
picture one considers only the lowest two stationary so- 
lutions of the original GP equation, (jTJ) , and self-trapping 
occurs due to a misalignment of their eigenenergies that 
suppresses the (Josephson) tunneling between the two 
states. The dynamics thus involves two quasi-degenerate 
states arising from the ground states of the unconnected 
wells. 



B. Beyond the usual two- mode dynamics 

What is new here is that near g\u — 2 (A ~ 60), the 
third eigenvalue becomes aligned with the second and, 
since the corresponding modes are located in different 
wells, tunneling of atoms is again allowed. But now be- 
tween modes 2 (left well) and 3 (right well). This cor- 
responds to the rise of A z in figure [T] beyond g\r> — 1. 
Figure 3a shows that by increasing the non-linear inter- 
action we go from the usual coupling between the two 
lowest modes, 1 and 2, see panel (c), to a coupling of the 
lowest mode of the left- well, 2, to the first mode of the 
right well, 3, see panel (d). This coupling, which is zero 
in absence of interaction, occurs due to the large nonlin- 
earity: it deforms the wave functions enough to enable 
the coupling between the ground state of one well and 
the first excitation in the other. It is a large and clearly 
density dependent effect which requires a change in the 
modes used and cannot be accommodated by varying the 
parameters of the usual two-mode picture. This effect is 
different in nature from what was reported in Refs. [l6r - 
[iH ] where the alignment takes place due to the presence 
of a large enough bias in the system. In our case, it is 
clearly a dynamical phenomenon which happens even for 
symmetric double-well potentials. The role played here 
by the non-linear interaction in modifying the single par- 
ticle states is more similar to the interaction blockade 
effect demonstrated for double-wells with few atoms in 
optical superlattices [l9j . 

Our result is also of different nature to the disappear- 
ance of self-trapping reported in Ref. [2(| ■ There, the 
authors explore the population of low-energy Bogoliubov 
excitations in the condensates of each of the wells finding, 
using a schematic model, a departure from self-trapping 
due to excitation of such low energy modes. By using 
the time-dependent GP equation to determine the con- 
densate wavefunction, the low energy Bogoliubov states 
are already incorporated into its time evolution. See Eq. 
(8.43) and the discussion in section VIII. E in Ref. [l|. 



Including them again along the lines of Ref. [20[ would 
be redundant. 

To further confirm our picture, in panel (b) of figure [3] 
we compare the frequencies of oscillation of z(t) 3 found 
in the full GP calculations, shown as empty circles, to the 
energy differences of the first three modes 4 . Up to g\u — 
1 the two mode description with states 1 and 2 works 
very well, but beyond that, the appropriate two mode 
model must involve states 2 and 3, and tunneling between 
wells increases again due to the progressive alignment 
between the energies of these two modes. The transition 
from (1,2) to (2,3) takes place smoothly, with the region 

1 9id 2 having more than one clear peak in the 
Fourier transform, z(lo), see panel (f) of Fig. [31 but an 
extremely small oscillation amplitude as seen in Fig. [1] 
In this region, the dynamics is governed by the first three 
modes coupled pairwise, (1,2) and (2,3). 

Above g\D — 2, the dynamics is dominated by modes 

2 and 3, depicted in panel (d) of Fig. [3] The strongest 
frequency extracted from the GP calculation now falls 
close to the E3 — E2 band (green). Following similar 
arguments to those in the derivation of Eqs. we can 
write down the new two mode equations: We denote the 
second and third modes of V c s(x, 0) as $l and A 
slight generalization of Eqs. ([2]) leads to 



i(i) = -2/C-v/l - z 2 (t) sin<ty(i) 



(6) 



6</>(t) = A + NUz{t) + 2K 



*(*) 



cos S(j)(t) 



where, 



A = (E° L - E° a ) + (U L - U R )/2 
U = (U L + U R )/2 



Ul [R) = A Jdx$i (R) (7) 
J dx ({l/2)d x ® L{R) d x $ L{R) + V(x 



E° L{R) = I dx 



Using gio = 1.2 to build the modes gives, A = —1.26, 
U = 0.35, and JC = —0.037. The prediction of this new 
two-mode model is shown by the dashed (blue) lines in 
Fig. [T] and Fig. [3Jb). As can be seen this (2,3) model 
works well in the range 1 < g\D < 3.5, giving a good ac- 
count of both the dominant frequency and the amplitude 
of the imbalance observed in the GP simulations. The 
transition from (1,2) to (2,3) coupling reflects also in the 
appearance of the node of the \ip(x,t)\ 2 near x = 2, ob- 
tained solving the GP equations as seen in Fig. |U This 



3 These frequencies correspond to those with the largest ampli- 
tudes in the direct Fourier transform |z(w)| of the function z(t), 
over one Rabi time. 

4 The energy of the fourth mode, see Fig. \3\ shows a constant in- 
crease, without avoided crossings in the range of gi£> considered, 
so that this mode remains mostly uncoupled to the rest and we 
will ignore it. 
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FIG. 4: Density of atoms along the x direction at a given time 
t — titabi/lO for several values of gin- The bands are built by 
joining the calculations at two different values of giD- 



is an observable feature which should be looked for ex- 
perimentally. 

Further increasing g±D, 9id ^ 6, our initial state has 
an average energy above the barrier, see dotted line in 
Fig. [3] (a), thus facilitating the flow of atoms between 
the wells. At high enough guj a certain equilibration 
of the imbalance can be expected, in line with Ref. [3] , 
mostly due to the sizeable contributions from modes with 
energies above the barrier. 



Critical g^ 3) for z(0) < 1 



Until now we have considered only cases with a max- 
imally imbalanced initial condition, z(0) = 1. For other 
initial conditions, z(0) < 1 we can estimate the critical 

(2 31 

value for the tunneling between modes (2,3), g\p , using 
similar arguments to the ones leading to Eq. Q. The 

value of g\p for z(0) = 1 corresponds to the one for 
which the second energy level (the first of the most pop- 
ulated well) equals the third energy level (the second of 
the less populated well). Thus fulfilling E\° + SE 2 = E 3 , 
which can easily be solved giving, g^^ — 2 3 / 4 y / 7r ~ 2.98, 
in good agreement with Figs. [2] and [3l 

Similarly, we can estimate the critical value, g^^\ for 
an initial state with a certain population imbalance zq = 
z(0) > 0. In this case, assuming the system remains 
mostly self-trapped, we can define the fraction of atoms 
on the left well, pl = (1 + zq)/2, and on the right well, 
Pr= (l — zo)/2. We can use a similar argument as before 
and assume the wave function of the system is essentially 
a harmonic oscillator wave function at each side of the 
trap. We assume, the left well stays mostly on the ground 
state, while the right side is promoted to the first excited 
one. Treating again the non-linearity as a perturbation, 



we have, (E x = E 1R , E 2 = E 1L , E 3 = E 2R , E 4 = E 2L ) 

5E 1 = giDPR(^ \\^ \ 2 \^°) = .giW(2 5/4 V^) , 
SE 2 = g 1D p L (rt°\\rt°\ 2 \rt°)=giDPL/(2 1/4 V^), 

SE 3 = giDP R (^ UZ \ 2 \i>2 ) = m)$Ei, 
SE^ = giDPLi^M ] 2 ]^ ) = (1/2)^2- 

The critical condition would correspond to have E 2 
that is E%° + SE 3 = E\° + SE 2 



_3_ 

71 



3 g\ c D PR. 

2 2 5 / 4 



1 

71 



2VVSF 



which gives, 



J2,3) 
'JlD 



4 2 3 / 4 



4p L - 3p R 



(8) 

(9) 
(10) 



,(2,3) 



2 3/4 



/5F[l + (7/4)p B ] 

,(2,3) 



For pr -C 1 we have, g 1D 

Thus, for z(0) close to 1, the critical value, g\ D °' , grows 
linearly with the initial fraction of atoms in the less pop- 
ulated well. 



D. Possible experimental conditions 

The phenomena presented in the previous sections rely 
on the adequateness of the Gross-Pitaevskii equation to 
describe the physics at large enough values of gio- The 
predicted reappearance of tunneling of atoms through 
the barrier beyond the usual self-trapped region as we 
increase g\o requires quantum fluctuations to be mini- 
mized during the whole process. Thus, we need to es- 
timate with realistic conditions whether the system can 
be brought to such gw values while remaining mostly 
condensed. 

We consider N atoms of 87 Rb trapped inside an axially 
symmetric trap, characterized by and its associated 
length a±_ = (hj 'A/cjjJ 1 / 2 . The double- well potential is 
built with a frequency uj x inside each well. The scaled 
ID coupling, q, can be written in the weakly interacting 
limit 5 as f 2 1.1 , 



<J 



4irH 2 a s 
2?ra 2 i M 



2ha s u± 



(11) 



To obtain meaningful estimates we consider the elon- 
gated trap conditions of Ref. [23[: uj x = 27r44.7 Hz and 
uj± = 27rll00 Hz. For these we can estimate the dilu- 
tion factor, rj = n!£>a s , where a s is the s-wave scattering 



As discussed in Ref. [2l|l , Eq. ((TJ can be obtained from the 3D GP 
one in the weakly interacting limit. For the strongly interacting 
limit the corresponding ID reduction would correspond to \4ff = 
V(x) + \o\/N\ip(x,t)\. We have checked that the phenomena 
described in the paper are also present in such reduction. 
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length and n±D is the one-dimensional density. In the 
Thomas- Fermi limit, this estimate is given by 



v = (am) 



2/.:! 



1 



Wl 3!/325/3 



(12) 



and for the trap conditions considered, r\ ~ 9 x 
10~ 3 (giD) 2 ^ 3 ■ Thus, an elongated trap ensures the di- 
luteness of the gas along the direction of the barrier, 
which is a necessary condition for the validity of the 
Gross-Pitaevskii approximation. The diluteness on the 
longitudinal direction also ensures the non-excitation of 
transverse modes of the elongated trap. This can be 
checked by comparing the transverse chemical potential 
with the transverse quantum, fku± . The transverse chem- 
ical potential can be written as, see Eq. (10) of [22j], 
fi± = hu±(l + 2rf) which is well below 2hu± for the val- 
ues of gio considered in this work provided the trap is 
cigar-shaped. 

This excitation of the third mode (the first excited 
state of the less-populated well) should show up macro- 
scopically in the experiments as an increase in the ampli- 
tude of oscillations of the population imbalance above a 
certain value of the interaction strength. Simultaneously, 
a node should appear in the center of the atomic cloud 
sitting in the less populated well, see Fig. 0J 

Experimentally, the challenge is twofold: first one 
needs to build a double-well potential with a single- 
particle spectrum close to the one described here, see 
Fig. Eh" This can be achieved by modifying the parame- 
ters of the optically produced double- well used in Ref. Q . 
Secondly, one needs to increase the number of atoms to 



have a non-linear term able to scan the transition de- 
scribed in this article. Such an experiment is within reach 
with current experimental setups [2^ |. 



IV. CONCLUSIONS 

We have demonstrated the possibility of exciting 
higher modes of the double-well potential through the 
dynamics. We have considered an initially imbalanced 
population and have shown that for a broad range of in- 
teraction energies the system remains self-trapped but 
not due to the dynamics of the two lower states of the 
Gross-Pitaevskii equation, as usually accepted, but of an- 
other mechanism involving a third state. This transition 
from the coupling between the first two (1,2) to the next 
two (2,3) states can be well characterized and understood 
by analyzing the static properties of the effective poten- 
tial (including interactions) which due to the fact that 
the system remains mostly self-trapped does not vary 
substantially with time. 
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